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INTRODUCTION 

Wave  scattering  from  a  sinusoidal  grating  has  been  a 
problem  of  considerable  interest  in  a  number  of  disciplines. 
It  is  of  interest  in  seismology  for  modeling  the  affects  of 
topography  on  seismic  wave  propagation.  Many  of  the 
discrete  wavenumber  modeling  techniques  that  have  been 
developed  for  sinusoidal  gratings  are  also  applicable  to 
interfaces  of  arbitrary  shape. 

In  this  paper  we  will  compare  results  and  numerical 
properties  of  four  different  discrete  wavenumber  modeling 
techniques.  All  of  the  comparisons  are  done  for  the  simple 
problem  of  a  sinusoidal  free  surface  grating  with  a 
vertically  traveling  plane  wave  incident  from  below  (Figure 
1)  .  Particle  motion  is  parallel  to  the  strike  of  the 
structure  (SH  waves) .  The  modeling  methods  we  will  examine 
are  the  Aki-Larner  method  (1970),  the  Wa':^  man  method 
(1975),  the  Waterman-Fourier  method  (Varaden,  et  al.  1987), 
and  the  Campillo-Bouchon  method  (1987).  We  will  refer  to 
these  methods  as  AL,  WR,  WF,  and  CB,  respectively.  All  of 
these  techniques  are  applicable  to  interfaces  of  arbitrary 
shape  and  can  in  principle  be  generalized,  to  fully  elastic. 


multi-layer,  and  3-D  models. 


A  solution  to  the  sinusoidal  grating  problem  was  1 


proposed  by  Payleigh  (1907) .  He  expanded  the  wavefield  in  - “ 

terms  of  up  and  down  going  plane  waves  and  solved  for  the  n / _ 
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condition  at  a  number  of  points  along  the  surface.  AL  is 
essentially  an  overdetermined,  least  squares  version  of  the 
Rayleigh  method,  formulated  in  the  wavenumber  domain. 
However,  since  Lippmann  (1953)  questioned  the  expansion  used 
in  the  Rayleigh  method,  the  so  called  "Rayleigh  ansatz" ,  the 
history  of  this  modeling  approach  has  been  shrouded  in 
controversy.  Lippmann  suggested  that  the  Rayleigh  expansion 
is  satisfactory  below  a  plane  bounding  the  topography,  but 
is  unsatisfactory  inside  the  the  peaks  (Figure  2) . 

Petit  and  Cadilhac  (1966)  proved  that  the  Rayleigh 
ansatz  could  not  be  valid  if  the  ratio  of  the  amplitude  to 
the  wavelength  of  the  sinusoid  (h/L)  is  greater  than  0.072. 
Millar  (1971)  proved  that  the  Rayleigh  ansatz  is  in  fact 
valid  if  h/L  is  less  than  0.072.  Other  authors  have  shown, 
e.g.  Wirgin  (1980)  and  Jiracek  (1973) ,  that  the  Rayleigh 
ansatz  is  usable  for  slopes  greater  than  the  Millar  limit. 
There  continues  to  be  heated  debate  over  this  issue.  (See 
Wirgin  1986,  and  Lakhtakia  et  al.,  1986.) 

In  addition  to  the  controversy  over  the  validity  of  the 
Rayleigh  ansatz ,  there  has  been  discussion  as  to  whether  the 
Waterman  method  (WR)  uses  the  ansatz.  Basically  this 
method,  also  known  as  the  T-matrix  method,  makes  use  of  the 
Helmholtz  formula  to  generalize  the  boundary  conditions. 
Waterman  originally  claimed  that  his  approach  was  not 
dependent  on  the  ansatz.  Today,  it  is  generally  accepted 
that  WR  uses  the  ansatz  and  is  therefore  limited  by  the 
slope  of  the  interface  (Lakhtakia  et  al.,  1985b). 
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WF  is  an  alternative  formulation  of  the  Waterman  method 
that  uses  the  expansion  proposed  by  Masel  et  al.  (1975) 
instead  of  the  Rayleigh  ansatz .  In  principle,  it  is  not 
limited  by  the  slope  of  the  interface. 

The  Campillo  and  Bouchon  method  (CB)  resembles  the 
aforementioned  techniques  in  some  ways,  and  it  is  clearly 
independent  of  the  Rayleigh  ansatz.  CB  is  a  collocation 
method  similar  to  AL,  but  it  parameterizes  the  reflected 
wavefield  in  terms  of  line  (or  point)  sources  distributed 
along  the  interface. 

We  will  compare  and  evaluate  the  four  techniques  as 
they  apply  to  the  sinusoidal  free  surface,  SH  wave  problem. 
In  particular  we  will  map  the  region  of  energy  conservation 
with  respect  to  slope  (h/L)  ,  frequency  (k*L),  and  N  (the 
number  of  plane  waves  used  in  the  expansion) .  We  will  also 
examine  the  wavefield,  at  the  interface  and  below  the  lower 
bounding  plane  of  the  topography. 
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METHODS 

In  this  section  a  unified  description  of  the  four 
methods  will  be  given  in  order  to  define  the  common 
features.  All  of  them  operate  on  time  harmonic  solutions  to 
the  wavefields,  and  the  time  dependence,  exp(iwt),  is 
suppressed.  In  all  four  methods  the  model  is  assumed  to  be 
periodic  in  the  x-direction,  leading  to  a  discrete  plane 
wave  representation  of  the  wavefield.  The  x  and  z 
wavenumber  components  are: 

kn  =  2nn/L  (1) 

and  7n  -  (w2/02  -  kn2)1/2  (2) 

where 

0  -  ( n/p)L/ i  -  shear  wave  velocity.  (3) 

Effects  of  adjacent  periods  can  be  suppressed  by  adding  a 
small  imaginary  component  to  w  which  can  be  defined  in  terms 
of  a  realistic  Q  parameter.  Each  of  the  methods  can  be 
stated  in  matrix  form. 

For  AL,  the  matrix  representation  is, 


where 


a  =  Q  r 


(4) 


rn  =  the  coefficients  in  the  Rayleigh  expansion  for 
the  reflected  waves,  i.e., 


Ur 


1  00 

X  rn 


i(knx  +  7nz) 
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z  >  £(x) 


(5) 


5 


am  =  the  Fourier  transform  of  the  traction  on  the 
interface  surface  due  to  the  incident  field. 
For  a  vertically  traveling  plane  wave  of  unit 
amplitude, 


a 


m 


fl 


1m7q/l 


i70?(x) 
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dx. 


(6) 


And 


fL 


mn 


-  -iM/L 


- 


a  i 

n*k 


km)x  irn€(x) 


dx . 


(7) 


Truncation  of  the  series  (5)  and  discretization  of  integrals 
(6)  and  (7)  leads  to  a  least  squares  solution, 

=  -M,  *  *  ‘  ,  0,  *  *  •  ,  +M  (8) 


where 


N  =  2M  +  1. 


(9) 


The  solution  to  (4)  minimizes  the  normal  stress  at  the 
surface  sample  points  in  a  least  squares  sense.  The  vector 
a  and  columns  of  the  matrix  G  can  be  calculated  using  the 
fast  Fourier  transform  ( FFT)  .  In  this  study  we  use  an  FFT 
of  length  128  and  solve  for  a  vector  r  of  length  N.  With  N 
<  128,  the  problem  is  overdetermined. 


The  Waterman  method  derives  from  the  Helmholtz  formula 
equations  (10  and  11),  aftet  plane  wave  expansions  have  been 
substituted  for  the  Green's  functions  and  the  wavefields. 
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Equation  10  represents  Huygen's  principle,  and  equation  11 
is  the  extended  boundary  condition.  (Waterman  1975) 


UR(z  >  O  ' 

► 

Uj(z  <  o  . 


l/4i 


*  -+■  oo 

n* [U+  Vq  -  g  VU+]  dx 

—  CO 


(10) 

(11) 


where 


and 


g  =  the  free  space  Green's  function. 


U+  =  the  wavefield  on  the  interface. 


Note  that  n-VU+  =  0  on  a  free  surface. 

Modeling  using  WR  or  WF  requires  that  we  solve  the 
following  coupled  set  of  matrix  equations. 


r  =  Q+  a  (12) 

b  =  Q-  a  (13) 


where 


and 


rn  = 


b  = 


the  coefficients  of  the  reflected  wavefield  as 
in  AL,  except  the  domain  is  limited  to  z  >  h, 


C*  = 


(14) 


the  coefficients  of  the  incident  wavefield,  i.e., 

1  n  -  0 

0  n  /  0 


the  coefficients  of  the  surface  field 
expansion . 


For  WR, 


If 


+M 
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ni-  -M 
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(15) 
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Q™  =  i/kL 


-  kmnx/nz>  e 


i(kmx  ±  7nn)  i(knx  -  7mO 

e  dx  (16) 


and  for  WF, 


U  =  I  Q 


+M  ik  x 

m 


m=-M 


m 


(17) 


Q__  =  i/kL 


"  kmnx/nz>  e 


i(V  1  iknx 


dx 


(18) 


Equations  12  and  13  correspond  to  equations  16  and  18, 
respectively,  with  an  implied  trucation  length  of  2M+1 
wavenumber  samples.  Equation  15  is  equivalent  to  the 
Rayleigh  ansatz  expansion,  albeit  expressed  in  terms  of 
upgoing  instead  of  downgoing  waves.  If  we  use  the 
alternative  surface  field  expansion  proposed  by  Masel  et  al. 
(1975) ,  the  matrix  equations  represent  the  Waterman  Fourier 
method  (WF) . 

It  is  interesting  to  note  that  although  they  are  used 
in  a  completely  different  context,  the  Q~  matrix  for  WF  is 
nearly  identical  to  the  transpose  of  the  G  matrix  used  in 
AL.  The  two  integral  expression  are  the  same  except  for  a 
factor  (nz)  in  the  integrand.  This  suggests  that  the 
elements  of  the  Q  matrix  can  be  generated  a  row  at  a  time 
using  an  FFT.  For  a  sinusoidal  interface,  the  integral  (18) 
simplifies  to  a  Bessel  function.  However,  in  this  study  we 
will  calculate  the  Q  matrices  using  FFT’s,  as  we  would  for 


an  interface  of  arbitrary  shape.  A  sample  length  of  128  is 
used,  the  same  as  with  AL. 

A  normalisation  scheme  can  be  applied  to  WR  or  WF  to 

improve  stability.  The  details  of  this  scheme  are  in 

Appendix  I . 

The  CB  method  is  a  collocation  method  like  AL.  Unlike 
the  other  methods,  however,  the  coefficient  vector  is  in  the 
space  domain  rather  than  the  horizontal  wavenumber  domain. 
The  wavefield  is  expanded  in  a  series  of  line  sources 
equally  distributed  in  x  along  the  interface: 

xn  -  (n-1)  L/N.  (19) 

For  CB  th*s  matrix  representation  is 

c  =  B  f  (20) 


where 


cn  =  the  traction  on  the  interface  due  to  the 
incident  field,  i.e.. 


1V<xn> 


c_  -  i  u  n_  7_  e 
n  ^  z  '  o 


(21) 


fm  =  the  line  source  coefficients,  i.e., 

IT!  '  ' 


+M  +M 


UR  =  l/(2i/*L)  l  f  l  (1/7D)  e 


iT'pl2  -  I  ikp(x  ‘  "Sn) 


m  ^  'p' 

m--M  p--M 


(22) 
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and 


+M  A 

Bnm  =  V(2iM)  £  (in-k/7  ) 
p=-M 


■m1 


ikp (xn  “ 


V 


(23) 


The  use  of  the  line  source  expansion  (22)  replaces  the 
Rayleigh  ansatz  and  Masel  expansions  used  in  'L  and  WF, 
respectively.  Equation  (20)  extinguishes  the  normal  stress 
at  the  source  points;  it  is  even-determined,  with  as  many 
equations  as  there  are  line  sources.  We  could  formulate  CB 
using  an  overdetermined  least  squares  solution  at  the  cost 
of  an  additional  loop  in  the  calculation  of  the  matrix 
elements.  This  addition  would  be  costly  in  terms  of 

execution  time. 
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RESULTS 

A  simple  way  to  evaluate  the  four  techniques  is  to  test 
for  energy  conservation.  In  each  numerical  experiment  we 
assume  that  a  monocroroatic,  vertically  traveling  plane  wave 
of  unit  amplitude  is  incident  on  a  sinusoidal  free  surface. 
Energy  conservation  requires  that  the  reflected  energy  flux 
equals  the  incident  energy  flux,  i.e.. 


+M 


l 

n«-M 


Re<V 


1, 


(24) 


as  used  by  others  (Larner  1970,  Lakhtakia  et  al.  1985a). 

The  interface  slope  is  described  by  h/L.  the  ratio  of 
the  amplitude  to  the  period  of  the  sinusoid.  The  frequency 
of  the  incoming  energy  is  described  by  k*L,  the  normalized 
frequency.  It  is  also  necessary  to  specify  N,  the  number  of 
coefficients  or  discrete  wavenumbers  samples  to  be  used  in  a 
particular  numerical  experiment.  The  region  of  energy 
conservation  in  the  two  dimensional  {N-(k'L) }  domain  has 
been  mapped  for  each  of  the  modeling  techniques  and  for 
various  interface  slopes. 

A  number  of  features  in  the  energy  conservation  maps 


are  worthy 

of 

mention.  As 

the 

normalized 

frequency 

increases , 

more 

coefficients 

are 

needed  to 

obtain 

convergence 

and 

energy  conservation. 

This  is 

expected 
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because  with  increasing  frequency  more  of  the  coefficients 
represent  propagating  waves.  The  WR,  WF,  and  AL  techniques 
are  all  limited  to  a  maximum  wavenumber  beyond  which  the 
solution  diverges  and  energy  is  not  conserved.  CB  does  not 
appear  to  have  such  a  limit. 

For  the  methods  that  are  wavenumber  limited,  the 
maximum  number  of  samples,  Nmax,  depends  on  frequency  and 
h/L,  the  slope  of  the  interface.  For  steep  slopes  (h/L  > 
0.15)  there  is  an  upper  limit  on  the  frequency  for  which 
convergence  occurs  at  any  wavenumber.  Therefore,  when  we 
discuss  slope  limitations  for  a  particular  modeling  method 
we  must  also  specify  frequency. 

Note  that  the  wavenumber  limit  is  a  numerical  feature 
not  necessarily  related  to  the  Rayleigh  ansatz .  For 
example,  the  WF  method  which  does  not  invoke  the  Rayleigh 
ansatz  has  an  NInax  for  slopes  both  above  and  below  the 
Millar  limit  of  h/L  =  0.072.  Also  the  the  WR  and  AL  methods 
have  wavenumber  limits  for  h/L  less  than  0.072. 

Further  scrutiny  of  Figure  3  reveals  that  Nmax  is 
always  in  the  evanescent  wave  region  where  coefficient 
amplitudes  increase  exponentially  with  z.  Use  of  these 
waves  in  the  solution  will  eventually  exceed  the  precision 
and  words ize  of  the  computer.  The  convergence  range  can  of 
course  be  extended  by  increasing  the  word  length  and 
precision. 

The  most  significant  feature  of  the  energy  conservation 
maps  is  that  the  WF  and  AL  methods  have  very  similar  regions 
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of  convergence.  For  steep  slopes,  the  two  methods  show 
convergence  up  to  about  the  same  maximum  frequency.  For  h/L 


=  0.1,  AL  and  WF  are  nearly  the  same  except  WF  diverges 
somewhat  more  rapidly  as  N  increases.  For  h/L  =  0.072, 
there  is  a  more  dramatic  difference  in  the  rate  of 
divergence.  Nonetheless,  these  data  support  the  observation 
made  by  Wirgin  (1980) ,  that  regardless  of  the  Rayleigh 
ansatz ,  the  WF  and  AL  methods  are  limited  to  the  same 
maximum  slope. 

The  performance  of  WR  falls  short  of  WF  or  AL.  For 
gradual  slopes,  WR  has  lower  Nmax  values,  and  for  steep 
slopes  WR  does  not  converge  at  all.  It  is  no  surprise  that 
WR  would  have  a  lower  Nmax  than  WF  since  Waterman's  original 
surface  field  expansion  introduces  an  additional  exponential 
factor  in  the  calculation  of  the  Q-matrix  elements. 
(Equation  15  as  compared  to  equation  17.) 

Use  of  the  normalization  technique  described  in 
Appendix  I  alters  the  region  of  convergence  for  the  WR  and 
WF  methods.  The  normalization  does  not  extend  the  maximum 
allowable  slope,  but  it  does  increase  Nmax,  allowing  more 
coefficients  to  be  used,  especially  at  frequencies  below  the 
maximum  for  a  particular  slope.  With  the  normalization,  the 
region  of  energy  conservation  for  the  WF  method  and  the  AL 
method  are  nearly  identical  even  at  low  slopes  (Figure  4)  . 
The  region  of  convergence  for  the  WR  method  is  also  extended 
by  use  of  the  normalization. 
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The  convergence  characteristics  of  the  CB  method  are 
substantially  different  from  the  other  three  methods.  For 
CB  there  is  apparently  no  limit  due  to  the  exponential 
factors  associated  with  evanescent  waves.  The  absolute 
value  in  the  formulation  of  the  CB  method  ensures  that  only 
exponentially  decaying  factors  are  in  the  matrix  elements. 
In  general,  more  wavenumber  samples  are  needed  with  the  CB 
method  to  obtain  the  same  accuracy  in  energy  conservation  as 
with  the  other  methods.  This  may  be  partially  due  to  the 
fact  that  CB  is  an  even-determined  collocation  method  as 
opposed  to  AL  which  is  over-determined. 

The  conservation  of  energy  constraint  is  one  way  to 
evaluate  the  four  modeling  techniques.  Another  such  test  is 
the  satisfaction  of  boundary  conditions.  Unfortunately,  for 
some  of  the  techniques  this  is  neither  practical  nor 
possible.  As  formulated,  the  CB  method  is  even-determined 
and  the  residuals  of  normal  stress  are  always  zero.  For  the 
Waterman  techniques  there  is  no  valid  expression  for  the 
stress  at  the  interface.  The  expansion  of  the  Green's 
function  in  terms  of  only  upgoing  or  downgoing  waves,  yields 
a  reflected  wavefield  solution  that  is  valid  only  in  the 
region  below  the  lower  bounding  plane. 

An  alternative  evaluation  for  the  modeling  techniques 
is  to  compare  the  displacements  at  the  interface.  It  is 
possible  to  calculate  displacements  at  the  interface  with 
all  of  the  methods,  even  the  Waterman  techniques  because 
equations  15  and  17  do  apply  immediately  at  the  surface. 
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Although  there  is  no  boundary  condition  for  displacements  to 
satisfy,  they  can  be  used  to  test  for  consistency  between 
the  methods.  Figure  5  shows  surface  displacements 
calculated  with  each  of  the  modeling  techniques  for 
interfaces  of  low,  moderate,  and  steep  slope.  In  each  case 
an  appropriate  frequency  and  number  of  samples  have  been 
chosen,  i.e.,  well  within  the  region  of  energy  conservation. 

For  small  slopes,  (h/L  =  0.1)  all  four  methods  agree 
well  in  terms  of  surface  displacement.  Note  that  the  AL  and 
WR  methods  agree  with  the  others  despite  the  fact  that  the 
slope  is  beyond  Millar's  theoretical  Rayleigh  limit  of 
0.072 . 

For  the  medium  and  steep  slope  models,  the  results  are 
different.  We  see  that  the  WF  and  CB  methods  agree  well  in 
terms  of  displacement  on  the  surface,  but  WR  and  AL  show 
marked  differences.  WR  does  not  satisfy  the  energy 
conservation  constraint  for  these  parameters,  so  it  is  not 
surprising  that  the  surface  displacements  are  different. 
For  AL  the  energy  conservation  constaint  is  satisfied,  but 
the  method  fails  to  yield  correct  surface  displacements  for 
slopes  significantly  beyond  the  Rayleigh  limit. 

It  is  useful  to  compare  the  wavefields  not  only  at  the 
surface,  but  also  below  the  lower  bounding  plane  of  the 
topography.  One  way  to  do  this  is  to  look  at  displacements 
at  the  lower  bounding  plane  (z  =  h)  .  Another  way  is  to 
compare  the  amplitude  spectra  of  the  reflected  energy. 
Figure  6  shows  the  displacements  at  z  =  h  for  the  same  slope 
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and  frequency  cases  for  which  we  examined  surface 
displacements.  For  h/L  =  0.1  all  four  methods  agree 
closely.  For  the  steeper  slope,  WF,  AL  and  CB  agree,  but  WR 
fails  to  give  a  consistent  result.  The  amplitude  spectra  in 
Figure  7  show  the  same  pattern  of  consistency.  For  steep 
slope  the  spectrum  for  WR  blows  up  in  the  evanescent 
wavenumber  range,  but  the  other  methods  agree  reasonably 
well.  (Note  the  logarithmic  scale.) 

For  slope  greater  than  the  Millar  limit  AL  fails  to 
yield  surface  displacements  consistent  with  the  other 
methods,  but  it  appears  to  give  the  correct  result  below  the 
lower  bounding  plane. 
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COMPUTATIONAL  CONSIDERATIONS 

In  general,  the  CB  method  requires  more  wavenumber  samples 
to  converge  and  is  slower  for  a  given  N,  than  the  other 
methods.  This  is  unfortunate,  since  it  performs  well  over  a 
broader  range  of  slope  and  frequency  parameters. 

Figure  8  shows  the  computation  time  for  each  of  the 
methods  on  the  CONVEX  Cl  supercomputer  using  single 
precision  arithmetic.  For  N  >  40  the  CPU  time  required  for 
CB  is  approximately  an  order  of  magnitude  greater  than  for 
AL  or  WF. 

Each  algorithm  consists  of  two  labor  intensive  parts, 
creating  the  N  X  N  matrix  and  solving  the  system.  The  AL, 
WR,  and  WF  algorithms  each  require  the  execution  of  three 
nested  loops  to  create  the  matrix.  Two  of  these  loops  can  be 
vectorized  by  the  CONVEX  fortran  compiler.  The  additional 
time  required  for  CB  is  due  to  an  additional  loop  implied  by 
equation  22  as  compared  to  equation  5  for  AL  and  equations 
15  and  17  for  WR  and  WF.  For  each  of  the  methods,  the 
LINPACK  Gaussian  elimination  routine  CGECO,  optimized  for 
the  CONVEX,  is  used  to  solve  the  system.  A  conjugate 
gradient  linear  equation  solver  was  tried,  but  proved  to  be 
slower  than  CGECO. 

In  Figure  8  we  see  the  dramatic  improvement  in  speed 
obtained  by  using  an  FFT  in  conjunction  with  the  Waterman 
Fourier  method. 
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8UMMARY 

We  have  presented  evidence  that  the  AL,  WR,  WF,  and  CB 
methods  are  useable  for  slopes  greater  than  the  theoretical 
limit  established  by  Millar  for  the  Rayleigh  ansatz.  For  a 
limited  class  of  problems,  (SH  waves  normally  incident  on  a 
sinusoidal  free  surface) ,  we  have  determined  the  maximum 
slope  and  frequency  that  each  method  can  handle.  This  has 
been  accomplished  by  mapping  the  region  in  slope/ frequency 
space  where  the  various  methods  converge  and  satisfy  the 
energy  conservation  requirement.  The  fact  that  the 
resulting  wavefields  are  consistent  between  methods 
indicates  that  the  convergent  solutions  are  indeed  valid. 

Of  the  four  methods  studied,  it  is  clear  that  CB  is  the 
most  stable  and  converges  over  the  broadest  range  of  slopes 
and  frequencies.  The  AL  and  WF  methods  have  convergence 
properties  very  similar  to  one  another  supporting  the 
conclusion  of  Wirgin  (1980)  that  neither  is  superior  for 
steep  slopes.  We  see  that  for  steep  slopes  ,  AL  gives 
erroneous  results  for  displacement  on  the  interface  but 
gives  valid  displacements  (consistent  with  WF  and  CB)  beyond 
the  topography.  WR  is  the  least  stable  of  the  four  methods 
and  converges  only  for  low  slopes. 

AL  possesses  a  potential  advantage  over  WF  for  multi¬ 
layer  problems  if  an  interface  is  within  the  topography  of 
an  adjacent  interface.  This  advantage  exists  only  for  low 
slopes  where  the  Rayleigh  expansion  is  valid  both  inside  and 
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outside  the  bounding  plane  of  the  topography,  i.e.  for  h/L  < 
0.072.  The  Waterman  methods  do  not  yield  solutions  to  the 
wavefield  inside  the  topography. 

We  have  introduced  a  normalization  technique  that 
improves  the  stability  of  WR  and  WF.  This  modification 
allows  more  wavenumber  samples  to  be  used  in  the 
calculations,  but  it  does  not  significantly  extend  the 
si ope/ frequency  limits  of  those  methods. 

In  terms  of  execution  time,  the  AL  and  WF  methods  are 
fastest.  Using  an  FFT  to  calculate  each  row  of  the  matrix 
in  WF  is  faster  than  computing  N  individual  integrals  or 
Bessel  functions.  For  nonsinusoidal  interfaces  the 
improvement  in  speed  will  be  most  dramatic. 

The  CB  method,  while  it  is  the  most  stable,  is  also  the 
slowest.  It  has  an  additional  loop  in  the  calculation  of 
the  matrix  elements  and  requires  more  wavenumber  samples  to 
achieve  convergence. 
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FIGURE  1  Schematic  of  the  numerical  experiments.  Uj  and  UR 
are  the  incident  and  reflected  wavefields.  Below  the 
interface  the  shear  modulus  and  density  are  p  and  p 
respectively. 


FIGURE  2  The  "Rayleigh  ansatz"  expands  the  reflected 
wavefield  below  the  interface  as  the  sum  of  downgoing  plane 
waves.  Lippmann  pointed  out  that  the  ansatz  was  incapable 
of  representing  upgoing  reflected  energy  in  the  region 
between  the  interface  and  the  plane  that  bounds  the 
topography . 


FIGURE  3  Energy  conservation  as  a  function  of  frequency  and 
number  of  coefficients.  Asterices  represent  reflected 
energy  flux  within  1%  of  incident  flux.  Plus  signs 
represent  1-5%  error.  Each  symbol  corresponds  to  a  separate 
numerical  experiment.  Each  frame  corresponds  to  a 
particular  technique  and  slope. 


FIGURE  4  Improvement  in  convergence  due  to  normalization. 
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FIGURE  5  Displacement  on  the  interface.  Each  frame 
corresponds  to  a  separate  numerical  experiment.  The  AL 
method  gives  anomalous  results  despite  satisfying  energy 
conservation. 


FIGURE  6  Dispacement  on  the  lower  plane  bounding  the 
topography  (z=h) ,  using  the  same  experimental  parameters  as 
in  Figure  5.  The  AL  method  agrees  with  WF  and  CB. 


FIGURE  7  Wavenumber  spectrum  of  reflected  energy  at  z=h, 
using  the  same  experimental  parameters  as  in  Figures  4  &  5. 


FIGURE  8  Computation  time  on  the  Convex  Cl  supercomputer  as 
a  function  of  N,  the  number  of  coefficients. 
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APPENDIX  I 


To  invoke  the  normalization  multiply  equations  12  and 
13  as  follow: 


W  r  =  W  Q+  a 
[W]"1  b  =  [W]-1  Q~  a 


where 


W 


mn 


-i7 

e 


n 


0 


n  -  m 
n  /  m 


Then  solve  the  new  system,  which  has  improved  stability: 


r '  =  Q+ '  a 

b'  =  Q~ 1  a 


r '  =  Q+,[Q",]”1b* 


Converting  to  the  original  basis  functions, 


r  =  W  Q+«  CQ"']"1  w  a. 
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Department  of  Geological  Sciences 
State  University  of  New  York 
at  Binghamton 
Vestal,  NY  13901 
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OTHERS  (United  States) 


Dr.  Monem  Abdel-Gawad 
Rockwell  International  Science  Center 
1049  Camino  Dos  Rios 
Thousand  Oaks,  CA  91360 


Prof.  Keiiti  Aki 

Center  for  Earth  Sciences 

University  of  Southern  California 

University  Park 

Los  Angeles,  CA  90089-0741 

Prof.  Shelton  S.  Alexander 
Geosciences  Department 
403  Deike  Building 
The  Pennsylvania  State  University 
University  Park,  PA  16802 

Dr.  Ralph  Archuleta 
Department  of  Geological  Sciences 
University  of  California  at  Santa  Barbara 
Santa  Barbara,  CA  93102 


Dr.  Thomas  C.  Bache,  Jr. 
Science  Applications  Inf  1  Corp. 
10210  Campus  Point  Drive 
San  Diego,  CA  92121  (2  copies) 


J.  Barker 

Department  of  Geological  Sciences 
State  University  of  New  York 
at  Binghamton 
Vestal,  NY  13901 

Dr.  T.J.  Bennett 
S-CUBED 

A  Division  of  Maxwell  Laboratory 
1 1800  Sunrise  Valley  Drive,  Suite  1212 
Reston,  VA  22091 

Mr.  William  J.  Best 
907  Westwood  Drive 
Vienna,  V A  22180 


Dr.  N.  Biswas 
Geophysical  Institute 
University  of  Alaska 
Fairbanks,  AK  99701 


Dr.  G.A.  Bollinger 
Department  of  Geological  Sciences 
Virginia  Polytechnical  Institute 
21044  Derring  Hall 
Blacksburg,  VA  24061 


Dr.  Stephen  Bratt 
Science  Applications  lnt’1  Corp. 
10210  Campus  Point  Drive 
San  Diego,  CA  92121 


Michael  Browne 
Teledyne  Geotech 
3401  Shiloh  Road 
Garland,  TX  75041 


Mr.  Roy  Burger 
1221  Serry  Road 
Schenectady,  NY  12309 


Dr.  Robert  Burridge 
Schlumberger-Doll  Research  Center 
Old  Quarry  Road 
Pidgefield,  CT  06877 


Dr.  Jerry  Carter 
Rondout  Associates 
P.O.  Box  224 
Stone  Ridge,  NY  12484 


Dr.  W.  Winston  Chan 
Teledyne  Geotech 
314  Montgomery  Street 
Alexandria,  V A  22314-1581 


Dr.  Theodore  Cherry 
Science  Horizon-.  Inc. 

710  Encinitas  Blvd.,  Suite  200 
Encinitas,  CA  92024  (2  copies) 


Prof.  Jon  F.  Claerbout 
Department  of  Geophysics 
Stanford  University 
Stanford,  CA  94305 


Prof.  Robert  W.  Clayton 

Seismological  Laboratory 

Division  of  Geological  &  Planetary  Sciences 

California  Institute  of  Technology 

Pasadena,  CA  91125 

Prof.  F.  A.  Dahlen 

Geological  and  Geophysical  Sciences 

Princeton  University 

Princeton,  NJ  08544-0636 


-4 


Prof.  Anton  W.  Dainty 
Earth  Resources  Lab 
Massachusetts  Institute  of  Technology 
42  Carleton  Street 
Cambridge,  MA  02142 

Prof.  Adam  Dziewonski 
Hoffman  Laboratory 
Harvard  University 
20  Oxford  St 
Cambridge,  MA  02138 

Prof.  John  Ebel 

Department  of  Geology  &  Geophysics 
Boston  College 
Chestnut  Hill,  MA  02167 


Eric  Fielding 
SNEE  Hall 
1NSTOC 

Cornell  University 
Ithaca,  NY  14853 

Prof.  Donald  Forsyth 
Department  of  Geological  Sciences 
Brown  University 
Providence,  RI  02912 


Prof.  An  Frankel 
Mai!  Stop  922 
Geological  Survey 
790  National  Center 
Reston,  VA  22092 

Dr.  Anthony  Gangi 
Texas  A&M  University 
Department  of  Geophysics 
College  Station,  TX  77843 


Dr.  Freeman  Gilbert 

Inst,  of  Geophysics  &  Planetary  Physics 
University  of  California,  San  Diego 
P.O.  Box  109 
La  Jolla,  CA  92037 

Mr.  Edward  Giller 
Pacific  Sierra  Research  Corp. 

1401  Wilson  Boulevard 
Arlington,  VA  22209 


Dr.  Jeffrey  W.  Given 
Sierra  Geophysics 
1  1255  Kirkland  Way 
Kirkland,  WA  98033 


Prof.  Stephen  Grand 
University  of  Texas  at  Austin 
Department  of  Geological  Sciences 
Austin,  TX  78713-7909 


Prof.  Roy  Greenfield 
Geosciences  Department 
403  Deike  Building 
The  Pennsylvania  State  University 
University  Park,  PA  16802 

Dan  N.  Hagedom 
Battelle 

Pacific  Northwest  Laboratories 
Battelle  Boulevard 
Richland,  WA  99352 

Kevin  Hutchenson 
Department  of  Earth  Sciences 
St.  Louis  University 
3507  Laclede 
St.  Louis,  MO  63103 

Prof.  Thomas  H.  Jordan 
Department  of  Earth,  Atmospheric 
and  Planetary  Sciences 
Massachusetts  Institute  of  Technology 
Cambridge,  MA  02139 

Robert  C.  Kemerait 
ENSCO,  Inc. 

445  Pineda  Court 
Melbourne,  FL  32940 


William  Kikendall 
Teledyne  Geotech 
3401  Shiloh  Road 
Garland,  TX  75041 


Prof.  Leon  Knopoff 

University  of  California 

Institute  of  Geophysics  &  Planetary  Physics 

Los  Angeles,  CA  90024 


Prof.  L.  Timothy  Long 
School  of  Geophysical  Sciences 
Georgia  Institute  of  Technology 
Atlanta,  GA  30332 


Prof.  Art  McGarr 
Mail  Stop  977 
Geological  Survey 
345  Middiefield  Rd. 
Menlo  Park,  CA  94025 
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Dr.  George  Mellman 
Sierra  Geophysics 
11255  Kirkland  Way 
Kirkland,  WA  98033 


Prof.  John  Nabelek 
College  of  Oceanography 
Oregon  State  University 
Corvallis,  OR  97331 


Prof.  Geza  Nagy 

University  of  California,  San  Diego 
Department  of  Ames,  M.S.  B-010 
La  Jolla,  CA  92093 


Prof.  Amos  Nur 
Department  of  Geophysics 
Stanford  University 
Stanford,  CA  94305 


Prof.  Jack  Oliver 
Department  of  Geology 
Cornell  University 
Ithaca,  NY  14850 


Prof.  Robert  Phinney 
Geological  &  Geophysical  Sciences 
Princeton  University 
Princeton,  NJ  08544-0636 


Dr.  Paul  Pomeroy 
Rondout  Associates 
P.O.  Box  224 
Stone  Ridge,  NY  12484 


Dr.  Jay  Pulli 
RADIX  System,  Inc. 

2  Taft  Court,  Suite  203 
Rockville,  MD  20850 


Dr.  Norton  Rimer 
S-CUBED 

A  Division  of  Maxwell  Laboratory 
P.O.  Box  1620 
La  Jolla,  CA  92038-1620 

Prof.  Larry  J.  Ruff 
Department  of  Geological  Sciences 
1006  C.C.  Little  Building 
University  of  Michigan 
Ann  Arbor,  MI  48109-1063 


Dr.  Richard  Sailor 
TASC  Inc. 

55  Walkers  Brook  Drive 
Reading,  MA  01867 


Thomas  J.  Sereno,  Jr. 

Science  Application  Int’l  Corp. 
10210  Campus  Point  Drive 
San  Diego,  CA  92121 


John  Sherwin 
Teledyne  Geotech 
3401  Shiloh  Road 
Garland,  TX  75041 


Prof.  Robert  Smith 
Department  of  Geophysics 
University  of  Utah 
1400  East  2nd  South 
Salt  Lake  City,  UT  84112 

Prof  S.  W.  Smith 
Geophysics  Program 
University  of  Washington 
Seattle,  WA  98195 


Dr.  Stewart  Smith 
IRIS  Inc. 

1616  North  Fort  Myer  Drive 
Suite  1440 

Arlington,  VA  22209 

Dr.  George  Sutton 
Rondout  Associates 
P.O.  Box  224 
Stone  Ridge,  NY  12484 


Prof.  L.  Sykes 

Lamom-Doherty  Geological  Observatory 
of  Columbia  University 
Palisades,  NY  10964 


Prof.  Pradeep  Talwani 
Department  of  Geological  Sciences 
University  of  South  Carolina 
Columbia,  SC  29208 


Prof  Ta-liang  Teng 
Center  for  Earth  Sciences 
University  of  Southern  California 
University  Park 
Los  Angeles,  CA  90089-0741 


Dr.  R.B.  Tittmann 

Rockwell  International  Science  Center 
1049  Camino  Dos  Rios 
P  O  Box  1085 
Thousand  Oaks,  CA  91360 

Dr.  Gregory  van  der  Vink 
IRIS,  Inc. 

1616  North  Fort  Myer  Drive 
Suite  1440 

Arlington,  VA  22209 

William  R.  Walter 
Seismological  Laboratory 
University  of  Nevada 
Reno,  NV  89557 


Dr.  Gregory  Wojcik 
Weidlinger  Associates 
4410  El  Camino  Real 
Suite  1 10 

Us  Altos,  CA  94022 

Prof.  John  H.  Woodhouse 
Hoffman  Laboratory 
Harvard  University 
20  Oxford  Street 
Cambridge,  MA  02138 

Dr.  Gregory  B.  Young 
ENSCO,  Inc. 

5400  Port  Royal  Road 
Springfield,  VA  22151-2388 
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GOVERNMENT 


Dr.  Ralph  Alewine  III 
DARPA/NMRO 
1 400  Wilson  Boulevard 
Arlington,  VA  01731-5000 


Mr.  James  C.  Battis 
GL/LWH 

Han  scorn  AFB,  MA  22209-2308 


Paul  Johnson 

ESS-4,  Mail  Stop  J979 

Los  Alamos  National  Laboratory 

Los  Alamos,  NM  87545 


Janet  Johnston 
GL/LWH 

Han  scorn  AFB,  MA  01731-5000 


Dr.  Robert  Blandford 
DARPA/NMRO 
1400  Wilson  Boulevard 
Arlington,  VA  87185 


Dr.  Katharine  Kadinsky-Cade 
GL/LWH 

Han  scorn  AFB,  MA  01731-5000 


Eric  Chad 
Division  9241 
Sandia  Laboratory 
Albuquerque,  NM  01731-5000 


Dr.  John  J.  Cipar 
GL/LWH 

Hanscom  AFB,  MA  01731-5000 


Mr.  Jeff  Duncan 
Office  of  Congressman  Markey 
2133  Rayburn  House  Bldg. 
Washington,  D.C.  20515 


Dr.  Jack  Evemden 
USGS  -  Earthquake  Studies 
345  Middlefield  Road 
Menlo  Park,  CA  94025 


Ms.  Ann  Kerr 
IGPP,  A-025 

Scripps  Institute  of  Oceanography 
University  of  California,  San  Diego 
La  Jolla,  CA  92093 

Dr.  Max  Koontz 
US  Dept  of  Energy/DP  5 
Forrestal  Building 
1000  Independence  Avenue 
Washington,  DC  20585 

Dr.  W.H.K.  Lee 

Office  of  Earthquakes,  Volcanoes, 
&  Engineering 
345  Middlefield  Road 
Menlo  Park,  CA  94025 

Dr.  William  Leith 
U.S.  Geological  Survey 
Mail  Stop  928 
Reston,  VA  22092 


Art  Frankel 
USGS 

922  National  Center 
Reston,  VA  22092 


Dr.  T.  Hanks 
USGS 

Nat'l  Earthquake  Research  Center 
345  Middlefield  Road 
Menlo  Park,  CA  94025 

Dr.  James  Hannon 

Lawrence  Livermore  Nat'l  Laboratory 
P.O.  Box  808 
Livermore,  CA  94550 


Dr.  Richard  Lewis 

Director,  Earthquake  Engineering  &  Geophysics 
U.S.  Army  Corps  of  Engineers 
Box  63 1 

Vicksburg,  MS  39180 

James  F.  Lewkowicz 
GL/LWH 

Hanscom  AFB,  MA  01731-5000 


Mr.  Alfred  Lieberman 
ACDA/Vl-OA'State  Department  Bldg 
Room  5726 
320 -21st  Street,  NW 
Washington,  DC  20451 
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Stephen  Mangino 
GL/LVVH 

Hanscom  AFB,  MA  01731-5000 


Dr.  Robert  Masse 
Box  25046,  Mail  Stop  967 
Denver  Federal  Center 
Denver,  CO  80225 


An  McGarr 

U.S.  Geological  Survey,  MS-977 
345  Middlefield  Road 
Menlo  Park,  CA  94025 


Richard  Morrow 
ACDA/VI,  Room  5741 
320  21st  Street  N.W 
Washington,  DC  20451 


Dr.  Keith  K.  Nakanishi 
LawTence  Livermore  National  Laboratory 
P.O.  Box  808,  L-205 
Livermore,  CA  94550 


Dr.  Carl  Newton 

Los  Alamos  National  Laboratory 

P.O.  Box  1663 

Mail  Stop  C335,  Group  ESS-3 
Los  Alamos,  NM  87545 

Dr.  Kenneth  H.  Olsen 

Los  Alamos  Scientific  Laboratory 

P.O.  Box  1663 

Mail  Stop  C335,  Group  ESS-3 
Los  Alamos,  NM  87545 

Howard  J.  Patton 

Lawrence  Livermore  National  Laboratory 
P.O.  Box  808,  L-205 
Livermore,  CA  94550 


Mr.  Chris  Paine 

Office  of  Senator  Kennedy,  SR  315 

United  States  Senate 
Washington,  DC  20510 

Colonel  Jeiry  J.  Perrizo 
AFOSR/NP,  Building  410 
Bolling  AFB 

Washington,  DC  20332-6448 


Dr.  Frank  F.  Pilotte 

HQAFTAC/TT 

Patrick  AFB,  FL  32925-6001 


Katie  Poley 
CIA-OSWR/NED 
Washington,  DC  20505 


Mr.  Jack  Rachlin 
U.S.  Geological  Survey 
Geology,  Rm  3  Cl 36 
Mail  Stop  928  National  Center 
Reston,  VA  22092 

Dr.  Robert  Reinke 
WL/NTESG 

Kirtland  AFB,  NM  87117-6008 


Dr.  Byron  Ristvet 

HQ  DNA,  Nevada  Operations  Office 

Attn:  NVCG 

P.O.  Box  98539 

Las  Vegas,  NV  89193 

Dr.  George  Rothe 
HQ  AFTAC/TGR 
Patrick  AFB,  FL  32925-6001 


Dr.  Michael  Shore 
Defense  Nuclear  Agency/SPSS 
6801  Telegraph  Road 
Alexandria,  VA  22310 


Donald  L.  Springer 

Lawrence  Livermore  National  Laboratory 
P.O.  Box  808,  L-205 
Livermore ,  CA  94550 


Dr.  Lawrence  Turnbull 
OSWR/NED 

Central  Intelligence  Agency,  Room  5G48 
Washington,  DC  20505 


Dr.  Thomas  Weaver 
Los  Alamos  National  Laboratory 
P.O.  Box  1663,  Mail  Stop  C335 
Los  Alamos,  NM  87545 


-9- 


J.J.  Zucca 

Lawrence  Livermore  National  Laboratory 
Box  808 

Livermore,  CA  94550 


GL/SULL 
Research  Library 

Hanscom  AFB  ,  MA  01731-5000  (2  copies) 


Secretary  of  the  Air  Force  (SAFRD) 
Washington,  DC  20330 


Office  of  the  Secretary  Defense 
DDR&E 

Washington,  DC  20330 


HQ  DNA 

Attn:  Technical  Library 
Washington,  DC  20305 


Defense  Technical  Information  Center 
Cameron  Station 

Alexandria,  VA  22314  (5  copies) 


Defense  Intelligence  Agency 
Directorate  for  Scientific  & 
Technical  Intelligence 
Washington,  DC  20301 


AFTAC/CA 

(STINFO) 

Patrick  AFB,  FL  32925-6001 


TACTEC 

Battelle  Memorial  Institute 
505  King  Avenue 

Columbus,  OH  43201  (Final  Report  Only) 


Mr.  Charles  L.  Taylor 
GULWH 

Hanscom  AFB,  MA  01731-5000 


DARPA/RMO/RETRIEVAL 
1400  Wilson  Boulevard 
Arlington,  VA  22209 


DARPA/RMO/Security  Office 
1400  Wilson  Boulevard 
Arlington,  VA  22209 


Geophysics  Laboratory 
Attn:  XO 

Hanscom  AFB,  MA  01731-5000 


Geophysics  Laboratory 
Attn:  LW 

Hanscom  AFB,  MA  01731-5000 


DARPA/PM 

1400  Wilson  Boulevard 

Arlington,  VA  22209 
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Dr.  Ramon  Cabre,  S.J. 
Observatorio  San  Calixto 
Casilla  5939 
La  Paz.  Bolivia 


Prof.  Hans-Peter  Harjes 
s  Institute  for  Geophysik 
Ruhr  University /Bochum 
P.O.  Box  102148 
x  4630  Bochum  1,  FRG 

Prof.  Eystein  Husebye 
NTNF/NORSAR 
P.O.  Box  51 

N-2007  Kjeller,  NORWAY 


Prof.  Brian  L.N.  Kennett 
Research  School  of  Earth  Sciences 
Institute  of  Advanced  Studies 
G.P.O.  Box  4 

Canberra  2601,  AUSTRALIA 

Dr.  Bernard  Massinon 

Societe  Radiomana 

27  rue  Claude  Bernard 

75005  Paris,  FRANCE  (2  Copies) 


Dr.  Pierre  Mecheler 
Societe  Radiomana 
27  rue  Claude  Bernard 
75005  Paris,  FRANCE 


Dr.  Svein  Mykkeltveit 
NTNF/NORSAR 
P.O.  Box  51 

N-2007  Kjeller,  NORWAY 
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CONTRACTORS  (Foreign) 
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FOREIGN  (Others ) 


Dr.  Peter  Basham 

Earth  Physics  Branch 

Geological  Survey  of  Canada 

1  Observatory  Crescent 

Ottawa,  Ontario,  CANADA  K1A  0Y3 

Dr.  Eduard  Berg 
Institute  of  Geophysics 
University  of  Hawaii 
Honolulu,  HI  96822 


Dr.  Michel  Bouchon 
I.R.I.G.M.-B.P.  68 
38402  St.  Martin  D'Heres 
Cedex,  FRANCE 


Dr.  Ililmar  Bungum 
NTNF/NORSAR 
P.O.  Box  51 

N-2007  Kjeller,  NORWAY 


Dr.  Michel  Campillo 
Observatoire  de  Grenoble 
I.R.I.G.M.-B.P.  53 
38041  Grenoble,  FRANCE 


Dr.  Kin  Yip  Chun 
Geophysics  Division 
Physics  Department 
University  of  Toronto 
Ontario,  CANADA  M5S  1A7 

Dr.  Alan  Douglas 
Ministry  of  Defense 
Blacknest,  Brimpton 

Reading  RG7-4RS,  UNITED  KINGDOM 


Dr.  Fekadu  Kebede 
Seismological  Section 
Box  12019 

S-750  Uppsala,  SWEDEN 


Dr.  Tormod  Kvaema 
NTNF/NORSAR 
P.O.  Box  51 

N-2007  Kjeller,  NORWAY 


Dr.  Peter  Marshal 
Procurement  Executive 
Ministry  of  Defense 
Blacknest,  Brimpton 

Reading  FG7-4RS,  UNITED  KINGDOM 

Prof.  Ari  Ben-Men  ahem 
Department  of  Applied  Mathematics 
Weizman  Institute  of  Science 
Rehovot,  ISRAEL  951729 


Dr.  Robert  North 

Geophysics  Division 

Geological  Survey  of  Canada 

1  Observatory  Orescent 

Ottawa,  Ontario,  CANADA  K1A0Y3 

Dr.  Frode  Ringdal 
NTNF/NORSAR 
P.O.  Box  51 

N-2007  Kjeller,  NORWAY 


Dr.  Jorg  Schlittenhardt 

Federal  Institute  for  Geosciences  St  Nat'l  Res. 
Postfach  510153 

D-3000  Hannover  51 ,  FEDERAL  REPUBLIC  OF 
GERMANY 


Dr.  Roger  Hansen 
NTNF/NORSAR 
P.O.  Box  51 

N-2007  Kjeller,  NORWAY 


Prof.  Daniel  Walker 
University  of  Hawaii 
Institute  of  Geophysics 
Honolulu,  HI  96822 


Dr.  Manfred  Henger 

Federal  Institute  for  Geosciences  8t  Nat'l  Res. 
Postfach  510153 
D-3000  Hanover  51,  FRG 


Ms.  Eva  Johannisson 
Senior  Research  Officer 
National  Defense  Research  Inst. 

P.O.  Box  27322 

S-102  54  Stockholm,  SWEDEN  -12- 


